Liquid–core cylindrical lens (LCL) is composed of two different cylindrical lenses stuck together as shown in Figs. 1(a) and 1(b), if two kinds of liquids with different refractive indexes (RIs) are injected into the liquid–core successively, the LCL will have the function of diffusion cell. When the collimated beams enter into the LCL along the path perpendicular to the axis of the lens, the image width at different heights has one-to-one sensitive dependence on the RI of the liquid at the corresponding height of the liquid core.^[1,2] The spatial RI resolution property of the LCL can be applied to the study of the liquid diffusion process,^[3,4] the liquid hygroscopic process,^[5] and the solid dissolution in a liquid.^[6] Diffusion coefficient is an important datum for studying liquid mass transfer, and its accurate measurement is of important scientific significance and application value for chemical engineering, biology, medicine, environmental protection, and other scientific fields.^[7,8] Diffusion coefficient is generally concentration-dependent. The traditional methods of measuring concentration-dependent diffusion coefficient require many diffusion experiments with different concentrations, which are time-consuming and tedious.^[9,10] In order to break these limitations, based on an LCL, a novel optical method of measuring rapidly and verifying reliably concentration-dependent liquid diffusion coefficient is demonstrated in this work. However, there are two challenges to be solved in measuring liquid diffusion coefficient by using the LCL. First, due to the diffusion process, the spatial profile of RI in the liquid core is not homogeneous and the ray’s propagation in the inhomogeneous medium does not follow the linear propagation property.^[11,12] Therefore, it is necessary to study the ray’s propagation and imaging law in the LCL composed of inhomogeneous media. Second, the liquid diffusion coefficient is dependent on the concentration in general, so that the diffusion equation has no analytic solution.^[13] The second problem to be solved in the work is how to use numerical calculation method to solve the diffusion equation, and obtain the concentration-dependent diffusion coefficient D(C) by comparing the calculated concentration profiles with experimental diffusion images.

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. (a) Three-dimensional (3D) imaging diagram of ray in homogeneous media, (b) 3D imaging diagram of ray in inhomogeneous media, (c) two-dimensional (2D) imaging diagram of ray in homogeneous media (xy plane), and (d) 2D imaging diagram of ray in inhomogeneous media (x’z plane).

In order to solve the first problem, the spatial and temporal RI profile n^s(z,t) of the diffusion solution is obtained by means of a special solution of the diffusion equation, i.e., an analytic solution as the diffusion coefficient is independent of concentration,^[1,13] together with the experimental relationship between solution concentration and its RI. Then, a nonlinear ray equation, which is a 2nd-order differential equation of ray propagation in inhomogeneous media, is established and solved numerically in liquid–core region^[14,15] Furthermore, when collimated beams pass through either the homogeneous (satisfying linear ray equation) or inhomogeneous solution (satisfying nonlinear ray equation),^[17–19] the ray tracing method is used to simulate the diffusion images on the focal plane of the LCL. The computational results of two kinds of diffusion images show that the imaging process in inhomogeneous solution can be treated approximately by the law of ray propagation in homogeneous solution, under the condition of small refractive index gradient of diffusion solution. Taking for example the diffusion process of triethylene glycol (TG) aqueous solution at room temperature (298 K), the conditions of small refractive index gradient are that the difference in initial concentration between two diffusion solutions is less than 60% and the diffusion time is greater than 270 min.

In order to solve the second problem, the experimental diffusion processes of TG aqueous solution under the initial concentration differences of 0%–60% and 40%–100% are studied in the LCL, and the corresponding experimental concentration spatial and temporal profiles C^e(z, t) of the diffusion solution are obtained during the diffusion time from 270 min to 540 min. Theoretically, the concentration-dependent diffusion coefficient is assumed to satisfy the polynomial expansion relation D(C) = D₀ × (1 + α₁C + α₂C² + α₃C³ + ⋅ ), where D is the diffusion coefficient in an infinite dilute solution which can be obtained by analyzing diffusion images in the range of lower concentration; α₁, α₂, and α₃ are under-determined parameters. Then, the finite difference method (FDM) is used to solve numerically the Fick diffusion equation under the initial and boundary conditions the same as the experimental conditions,^[20,21] the calculated concentration profiles (Cⁿ(z,t)) by varying the under-determined parameters are compared with the profile C^e(z,t), the parameters (α₁, α₂, α₃, …) corresponding to the minimum concentration standard deviation between Cⁿ(z,t)s and C^e(z,t) are selected to determine D(C).

To verify the correctness of the obtained D(C), on the one hand, the D(C) is compared with the literature values measured by interferometry at the same temperature.^[22,23] On the other hand, the obtained D(C) is substituted into the diffusion equation, the numerical solutions of spatial and temporal RI profiles nⁿ(z_j, t_i) are calculated by means of the FDM. Based on the ray propagation law in inhomogeneous media, the ray tracing method is again used to calculate the dynamic images of whole experimental diffusion process (from 25 min to 540 min). The simulated results are in good agreement with the experimental images in image contour, focal position, and light intensity distribution, which verifies not only the correctness of D(C) measurement results, but also the reliability of treatment on ray propagation in the LCL composed of inhomogeneous diffusion solution in this work.

If two kinds of liquids with different values of RI of n₄ and n₁ (n₄ > n₁, the height of each liquid column is H) are injected into the liquid core successively, and a complementary metal–oxide–semiconductor (CMOS) detector is placed at the position where the RI (n_c = n₁) of upper liquid can be sharply imaged, thus, the liquid below presents a diffuse image as shown in Fig. 1(a). As the two liquids diffuse to each other, the solution forms a dynamic RI gradient distribution (n₄ > n₃ > n_c > n₂ > n₁) along the diffusion direction (z axis). If the detector is at the focal position of the LCL (f = f_c, n = n_c), according to the linear propagation law of ray in homogeneous media, the collimated beams entering into the LCL on the xy plane keep propagating in the same plane as shown by the dotted line in Fig. 1(b). The diffusion image is sharply imaged at the position z = z_c, while the detector presents a “beam waist” shaped dispersion image at other locations of z ≠ z_c. As shown on the right of Fig. 1(b), its image width (W_j) at the height of z_j has a one-to-one sensitive dependence on the RI (n = n_j) of the solution at the corresponding height z_j, which can be uniquely determined by the geometric relation and imaging law shown in Fig. 1(c), or accurately determined by experimental method.^[24]

However, according to the ray propagation law in inhomogeneous media, collimated beams passing through the liquid core will deviate from the entrance plane (xy plane) and propagate as shown by the solid line in Fig. 1(b) due to the fact that the RI of diffusion solution in liquid core is inhomogeneous along the z axis. Can the imaging law in homogeneous media be used to treat approximately the imaging process in inhomogeneous media, and can the light beam deviation from entrance plane be ignored, under the condition that the spatial and temporal RI profile of diffusion solution n(z, t) is undetermined? This is what needs to be analyzed next.

Based on the Snell law, the ray tracing method is used to track the position of each collimated ray after passing through the LCL and then hitting on a virtual CMOS detector. Light intensity profile is obtained by counting up the ray numbers reaching each pixel of the detector. As shown in Fig. 1(c), the vertices of the four spherical surfaces of the LCL intersect the x axis at a₁, a₂, a₃, and a₄, and their curvature radii are R₁, R₂, R₃, and R₄, respectively. The lens material is K9 glass (n = 1.5163). Taking for example one ray incident at point P₁(x₁, y₁, z₁), its ray equations in different light paths are established in the following.

After entering into the LCL through the 1st spherical surface along the xy plane, the ray keeps propagating on the xy plane (homogeneous media, n = n₀) as shown in Figs. 1(c) and 1(d). The coordinates of P₁ can be calculated by the spherical surface equation and the ray equation of the incident ray y = –h. The slope of the normal line at P₁ is , the incident and refractive angles on the 1st spherical surface are θ₁ = arcsin (h / R₁ ) and , respectively. The refractive ray intersects with the 2nd spherical surface at P₂ , and its slope is . Therefore, the coordinates of P₂ can be calculated by the linear ray equation y = k₁(x – x₁) + y₁ and the 2nd spherical surface equation .

Due to an RI inhomogeneous distribution along the z axis in the liquid core, when the ray P₁P₂ enters into the liquid core through the 2nd spherical surface, refractive ray P₂P₃ will deviate from the xy plane propagate on the plane (inhomogeneous media), and intersects with the 3rd spherical surface at P₃, as shown in Figs. 1(c) and 1(d).

The trajectory component of the refractive ray P₂P₃ in the xy plane is , the slope of the normal line at is , the incident and refractive angles on the 2nd spherical surface are

To express the ray equation on the plane, the dynamic coordinate system x’z is established as shown in Fig. 1(d). In the x’z coordinate system, the coordinates of P₂ are and z = z₂, and the line determines the x’ axis. The path infinitesimal element of a ray propagating in the x’z plane satisfies the geometric relationship as shown in Fig. 1(d)

Equation (3) has no analytic solution in general, but its numerical solution can be calculated on condition that the RI function n(z, t) is known. For example, if the term on the right-hand side of Eq. (3) is equal to the value of 2u at z = z₁ and t = t₀, the nonlinear ray equation gives out a quadratic curve on the x’z plane, which satisfies

As the ray leaves the liquid core from the 3rd spherical surface, refractive ray P₃P₄ will propagate on the plane (homogeneous media, n = n), and intersect with 4 th spherical surface at P₄, as shown in Figs. 1(c) and 1(d). The trajectory component of P₃P₄ in the xy plane is , the slope of the normal line at is , the incident and refractive angles on the 3rd spherical surface are

To express the ray equation on the plane, the dynamic coordinate system x’z is established on the plane as shown in Fig. 1(d). In the x’z coordinate system, the coordinates of P₃ are and z = z₃, and the line determines the x’ axis. The ray equation of P₃P₄ is , where . Substituting into z (x’), . The refractive angle can be calculated by the incident angle and refraction law.

3.4. Linear path of P₄P₅

After the ray leaves the 4th refractive spherical surface, refractive ray P₄P₅ will propagate on the plane (homogeneous media, n = 1), and intersect with CMOS detection at P₅ as shown in Figs. 1(c) and 1(d). The trajectory component of the refractive ray P₄P₅ in the xy plane is . The slope of the normal line at is . The incident and refractive angles on the 4th spherical surface are and , respectively. The satisfies the linear ray equation y = k₄(x – x₄) + y₄, where the slope k₄ of is

Therefore, the coordinates of , y₅) can be calculated by the linear ray equation and CMOS plane x = x₅.

Like P₃P₄, the linear ray equation of P₄P₅ in x’z coordinate system is , where , is given by . Hence, the coordinates of P₅(x₅, y₅, z₅) on the CMOS plane are finally determined by the ray tracing method from the ray at incident point P₁(x₁, y₁, z₁).

Based on the analysis in Section 4, under the conditions of initial concentration difference C₂ – C₁ = 0.6 and room temperature (298 K), two groups of diffusion experiments of TG aqueous solution are carried out to combine a complete concentration range from C = 0 to 1, the C₁ = 0, C₂ = 0.6 (0%–60%) and C₁ = 0.4, C₂ = 1.0 (40%–100%) for the first and second groups, respectively. Taking the first group for example, a digital injection pump is used to inject 60% TG aqueous solution (C₂ = 0.6) into the lower half of the LCL (H = 25 mm). Waiting for 10 min to reduce liquid turbulence interference, distilled water (C₁ = 0) is slowly injected into the upper half of the LCL at a speed of 2.0 ml/min (H = 25 mm), ensuring that there is no significant convection interference between the two liquids. The time when the two liquids contact each other is set to be the initial diffusion time (t = 0). The CMOS detector is fixed at the position where the collimated beam can be sharply imaged in the liquid thin layer with RI (n) = n_c = 1.3482, and images are collected in the interval of 2 min. The experimental operation for the second experimental group is similar to that for the first group, but the CMOS detector is fixed at the position where the collimated beams could be sharply imaged in the liquid thin layer with RI (n) = n_c = 1.3994.

To simultaneously satisfy the diffusion time required by the above simulated calculation and the initial conditions in Eq. (6), diffusion images during 306 min–380 min and 450 min–540 min are selected to calculate the concentration-dependent diffusion coefficients D(C) for the first and second diffusion group, respectively.

The typical diffusion images are shown in Fig. 3, where figures 3(a) and 3(b) show the diffusion images of the first group at t₀ = 320 min and 380 min. Figures 3(c) and 3(d) show the diffusion images of the second group at t₀ = 450 min and 540 min. After the processing of the binarization digital image, the extraction of the image width, and the transformation of the image width and the solution concentration, a series of experimental concentration spatial and temporal profiles C^e(z_j, t)s are achieved. Some of curves corresponding to those in Fig. 3 are plotted by the rectangular points as shown in Fig. 4.

Fig. 3.

	Figure Option ViewDownloadNew Window
	Fig. 3. Typical experimental diffusion images for [(a) and (b)] 0%–60% diffusion group at t = 320 min (a) and 380 min (b), and [(c) and (d)] 40%–10% diffusion group at t = 450 min (c) and 540 min (d).

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. Spatial and temporal concentration profiles Ce(zj, t)s corresponding to Fig. 3, at t = 320 min (a), 380 min (b), 450 min (c), and 540 min (d), with rectangular points representing experimental values, and solid red lines denoting calculated values.

In general, the diffusion coefficient is dependent on concentration and written as D(C). The D(C) can be expressed in the form of polynomial as

Partial differential equation (9) has no analytic solution on condition that the D(C) is undetermined. The algorithm of FDM is used to calculate the numerical solution of Eq. (9) in this paper. The spatial distance is divided into discrete M + 2 terms, j = 0 and M + 1 are the boundary terms. Equations (9) and (6) are changed into the difference quotient form, respectively,

Equation (10) is further expanded into M linear equations along the diffusion direction (z = jh, j = 0,1,2,…, M + 1). By varying under-determined parameters [(α₁)_k, (α₂)_k, (α₃)_k], a series of calculating Cⁿ(z_j, t_i)s of Eq. (10) at a special moment t_i = t is used to compare with the experimental profile C^e(z_j, t₀), and standard deviation σ_k is calculated. The minimum value of σ_k is the best-fit parameters of TG at t_i = t₀, the corresponding concentration values calculated are shown in Fig. 4 by the solid lines. and C^e(z_j, t₀) are in good agreement in the whole concentration area, which indicates that the D(C) relation calculated by comparing the experimental curve is reliable.

For the first diffusion experimental group, 10 diffusion images acquired in the interval of 5 min during the period from 305 min to 350 min, are used to calculate D(C) relations, the average values are listed in Table 1, which satisfy D₁(C) = 0.730 × 10⁻⁵(1 – 0.916C – 0.058C² – 0.045C³) cm²/s in the concentration range from C = 0.0 to 0.6, and D₁(C) versus mass fraction C is shown in Fig. 5 by the red empty circles. For the second group, 10 diffusion images acquired in the interval of 10 min during the period from 450 min to 540 min, are used to calculate D(C) relations, the average values are listed in Table 2, which satisfy D₂(C) = 0.730 × 10⁻⁵(1 – 0.860C – 0.008C² – 0.006C³) cm²/s in the concentration range from C = 0.4 to 1.0, and the D₂(C) versus mass fraction C is shown in Fig. 5 by the red solid circles.

Fig. 5.

Figure Option ViewDownloadNew Window

Fig. 5. Comparison of D(C) between measured values and literature values. Red empty circle represents measured value of D1(C) in concentration range of 0%–60%. Solid red circle denotes measured value of D2(C) in concentration range of 40%–100%. Solid line refers to fitting curve D(C) in concentration range of 0%–100%. Upper triangle represents measured values of literature 1. Lower triangles are for measured values of literature 2. Error bars, which are magnified 5 times for the sake of clarity, are determined by ten-time independent measurements.

Table 1.

Table 1.

Table 1. Measurement and calculation results of D(C) in a range of 0%–60% concentration. . D(C)/10−5 cm2⋅s−1 C/MF 0.0 0.1 0.2 0.3 0.4 0.5 0.6 α1 α2 α3 0.730 0.663 0.594 0.525 0.454 0.381 0.307 –0.916 –0.058 –0.045 STD 0.00 0.001 0.002 0.003 0.005 0.007 0.010 0.008 0.019 0.022 RSD/% 0.00 0.151 0.337 0.571 1.101 1.837 3.257 D1(C) = 0.730 × 10−5(1 – 0.916C – 0.058C2 – 0.045C3) cm2/s

Table 1. Measurement and calculation results of D(C) in a range of 0%–60% concentration. .

Table 2.

Table 2.

Table 2. Measurement and calculation results of D(C) in a range of 40%–100% concentration. . D(C)/10−5 cm2⋅s−1 C/MF 0.4 0.5 0.6 0.7 0.8 0.9 1.0 α1 α2 α3 0.478 0.414 0.350 0.286 0.222 0.157 0.092 –0.860 –0.008 –0.006 STD 0.002 0.003 0.004 0.004 0.004 0.005 0.006 0.006 0.003 0.003 RSD/% 0.418 0.725 1.143 1.399 1.802 3.185 6.522 D2(C) = 0.730 × 10−5(1 – 0.860C – 0.008C2 – 0.006C3) cm2/s

Table 2. Measurement and calculation results of D(C) in a range of 40%–100% concentration. .

To simply express the liquid diffusion coefficients D(C) in the whole concentration range by one formula, let the diffusion coefficients be D₁(C) in the range of C = 0 to 0.4, D₂(C) in the range of C = 0.6 to 1.0, and average value of D₁(C) and D₂(C) in the range of C = 0.4 to 0.6, the D(C) relation is refitted by the least square method in the range of C = 0 to 1, the fitted result is

It is worth noting that equation (12) is reliable in the whole diffusion process, though it is obtained by using the images acquired during the diffusion time from 305 min to 540 min. A prominent advantage of measuring D(C) relation introduced in the paper is that the correctness of deduced D(C) can be verified by comparing simulated images with experimental images in the whole diffusion process.

The algorithm of FDM is used to directly solve the diffusion Eq. (5) under the conditions of Eq. (11) and deduced D(C), and obtain the numerical solution of spatial and temporal concentration profiles Cⁿ(z_jt_i) (j = 0,1,…, M + 1; i = 0,1,2,…). The Cⁿ(z_jt_i) are then converted into spatial and temporal RI profiles n(z_jt_i) (j = 0,1,…, M + 1; i = 0,1,2,…) by means of experimental relationship of n(C) = 0.1251C + 1.3328. The undetermined coefficient u in nonlinear optical Eq. (4) satisfies the difference form

Fig. 6.

	Figure Option ViewDownloadNew Window
	Fig. 6. Comparison between experimental diffusion images and simulated images at RI = nc = 1.3482 for the first diffusion group (0%–60%): [(a′), (b′), (c′)] experimental diffusion images at t = 25, 210, and 380 min; [(a), (b), (c)] simulated images in “homogeneous solution”; [(a″, b″, c″)]: simulated images in inhomogeneous solution. Longitudinal scales represent solution concentrations.

Fig. 7.

	Figure Option ViewDownloadNew Window
	Fig. 7. Comparison between experimental diffusion images and simulated images at RI = nc = 1.3994 for the second diffusion group (40%–100%): [(a′), (b′), (c′)] experimental diffusion images at t = 40, 360, and 540 min; [(a), (b), (c)] simulated images in “homogeneous solution”; [(a″), (b″), (c″)] simulated images in inhomogeneous solution. Longitudinal scales represent solution concentrations.

When the virtual CMOS detector is placed at the position where the thin layer of RI = n = 1.3482 can be sharply imaged, the experimental diffusion images for the first group at t = t₀ = 25, 210, and 380 min are shown in Figs. 6(a′)–6(c′). The two kinds of simulated diffusion images are also shown in Fig. 6 for the convenience of comparison, where figures 6(a)–6(c) are the images ignoring the ray deviation from entrance plane and figures 6(a″)–6(c″) are the images that the ray deviation is calculated strictly when light ray passes through the inhomogeneous solution in the used LCL.

The comparison of the experimental diffusion images with two kinds of simulated diffusion images clearly shows some points as follows. (i) Under strict simulation conditions, the simulated diffusion images shown in Figs. 6(a″)–6(c″) are in very good agreement with the experimental images shown in Figs. 6(a′)–6(c′) in a wide diffusion time range from t = 25 min to 380 min, whether on the image contours, focal positions, or intensity distributions. (ii) Ignoring the ray deviation from entrance plane, only at the later stage of the diffusion process (t = 380 min), does the simulated diffusion image figure 6(c) coincide with the experimental diffusion image (Fig. 6(c′)).

As the virtual CMOS detector is placed at the position where the thin layer of RI = n = 1.3994 can be sharply imaged, the experimental diffusion images for the second group at t = t = 40, 360, and 540 min are shown in Figs. 7(a′)–7(c′). The two kinds of simulated diffusion images are also shown in Fig. 7 for the convenience of comparison, where figures 7(a)–7(c) show the images ignoring the ray deviation from entrance plane. Figures 7(a″)–7(c″) show the images that the ray deviation is calculated strictly when ray passes through the inhomogeneous solution in the used LCL.

The same phenomena appearing in Fig. 6 are also found in Fig. 7, implying that the location of virtual CMOS detector does not affect the simulation calculation of diffusion image. In addition, it is interesting to note that the focal positions of the experimental diffusion images, marked by the dotted line arrows both in Figs. 6 and 7 drift slowly upward with the diffusion time going by. This macroscopic image movement is caused by microscopic molecular thermal movement, and the simulated diffusion images reflect faithfully the law of thermal movement.

Furthermore, as the location of virtual CMOS detector moves to the position where the layer of n_c = 1.3398 is sharply imaged, the whole dynamic diffusion process is simulated and compared with the experimental diffusion process. The comparison between two diffusion processes is shown in visualization 1, which not only provides a vivid observation of the dynamic diffusion process, but also shows a solid verification of deduced D(C) due to the high consistency between the experimental process and the simulated process.

Ray equations in an LCL are established by using the ray tracing method, and the diffusion images on the focal plane of the LCL are simulated numerically, showing that the ray propagation law in homogeneous media can be used to calculate approximately the diffusion images in inhomogeneous media under the condition of small RI gradient of diffusion the solution. Based on the calculated conditions, the liquid diffusion process of TG aqueous solution in the LCL is experimentally implemented at room temperature, thus obtaining the liquid diffusion coefficient of TG aqueous solution D(C), which is dependent on the solution concentration. The spatial and temporal concentration profile C^e(z,t) is acquired by selecting the experimental diffusion images at the late stage of the diffusion process. The D(C) is expressed in the form of polynomial with a series of underdetermined parameters. The FDM method is used to calculate the spatial and temporal concentration profiles Cⁿ(z,t) of diffusion solution by varying the underdetermined parameters, and D(C) = 0.730 × 10⁻⁵(1 – 1.068C + 0.444C² – 0.2465C³) cm²/s is obtained by comparing the C^e(z,t) and Cⁿ(z,t). Based on the ray equations and the RI spatial and temporal profiles n(z, t) in the whole diffusion process calculated by the obtained D(C), we simulate strictly the diffusion images in an inhomogeneous solution. The simulated results are highly consistent with the experimental images in image contour, focal position, and light intensity distribution, which verifies not only the correctness of D(C) measurement result, but also the reliability of studying ray propagation in the LCL composed of inhomogeneous diffusion solution.